Révéler les pulsations urbaines circadiennes à l’aide d’un cartogramme animé : l’exemple de Bogotá

Définition de l’environnement de travail

Chargement des packages

    library(downloader) # pour télécharger des fichiers depuis un serveur
    library(tidyverse) # pour simplifier la syntaxe
    library(sf) # pour manipuler des objets SimpleFeature
    library(sp) # pour manipuler des objets Spatial*DataFrame
    library(spatstat) # pour calculer le lissage spatial
    library(maptools) # pour convertir l'objet sp en objet ppp (spatstat)
    library(cartography) # pour cartographier les résultats
    library(cartogramR) # pour créer les cartogrammes
    library(raster)  # pour manipuler les rasters
    library(terra) # pour manipuler les rasters
    library(stars) #pour la vectorisation du raster
    library(rmapshaper) # pour généraliser les contours (simplification des géométries)
    library(dplyr) #pour manipuler des tableaux de données
    library(openxlsx) #pour lire des fichiers EXCEL
    library(ggplot2) #pour tracer des graphiques esthétiques
    library(tidyr) #pour manipuler des matrices
    library(rgdal) # pour manipuler des objets ayant une composante spatiale
    library(rgeos) # pour calculer des centroides
    library(spdep) # pour calculer l'auto-corrélation spatiale
    library(geoR) # pour calculer le semi-variogramme empirique
    library(magick) # pour créer un gif animé
    library(av) # pour créer une vidéo mp4
    library(knitr) # pour afficher les images exportées
    library(filesstrings) # pour déplacer un fichier d'un dossier à un autre

Téléchargement et décompression du jeu de données

    download("https://github.com/ESO-Rennes/Animated-Cartograms/raw/main/Data_Pulsation.zip", destfile=paste0(getwd(),"/Data_Pulsation.zip"), mode="wb", overwrite = TRUE)
    unzip(zipfile = "Data_Pulsation.zip", exdir = ".")

Chargement de la couche SIG (fond de carte)

    EMU2019 <- st_read(dsn = "Data_Pulsation", layer = "EMU2019", stringsAsFactors = FALSE)

## Reading layer `EMU2019' from data source 
##   `F:\Rennes\Recherche\ANR MoDurAL\Valorizacion\IJC\version anglaise\Derniers_traitements\Envoi Florent 12 01\Data_Pulsation' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 56 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N

Chargement du tableau contenant les données de l’enquête

    viajes <- read.xlsx("Data_Pulsation/ViajesEODH2019.xlsx") #134497 observations

Définition d’une fonction pour rendre les couleurs transparentes (utile pour le fondu enchaîné)

    t_col <- function(color, percent = 50, name = NULL) {
      rgb.val <- col2rgb(color)
      t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
                   max = 255,
                   alpha = (100 - percent) * 255 / 100,
                   names = name)
    }  
    t_col_palette <- function(palette, percent = 50) {
      palette_trans <- palette
      for (j in 1:length(palette_trans)){
        palette_trans[j] <- t_col(palette_trans[j], percent = trans)
      }
      return(palette_trans)
    }

Création de la table de stock de population par UTAM avec pas de temps de 15 minutes

Calcul des durées de trajet (prise en compte des trajets qui commencent un jour et terminent le suivant)

    viajes$duracion <- viajes$p31_hora_llegada - viajes$hora_inicio_viaje

    for(i in 1:nrow(viajes)){
        if(viajes$p31_hora_llegada[i]<viajes$hora_inicio_viaje[i]){
        viajes$duracion[i] <- viajes$duracion[i]+1
        }
    }

Calcul du temps de trajet moyen en heures. Le temps de trajet moyen (aller-simple) est de 50 minutes

    avg = 24*mean(viajes$duracion) 

Tri pour écarter les trajets aberrants (0 minute)

    viajes2 <- viajes %>% filter(duracion>0)

Traitement des trajets dont le point de départ ou d’arrivée se trouve en dehors de Bogotá DC. On crée une UTAM virtuelle regroupant les zones hors DC.

    #Repérage des trajets qui sortent ou entrent dans le DC 
    excl <- c("N/A", "UPR1", "UPR2", "UPR3")

    viajes2$utam_destino[is.na(viajes2$utam_destino)] <- "N/A"
    viajes2$utam_origen[is.na(viajes2$utam_origen)] <- "N/A"

    # création d'un UTAM virtuel "UTAM99"
    for (i in 1:nrow(viajes2)){
      if (viajes2$utam_destino[i] %in% excl){
        viajes2$utam_destino[i] <- "UTAM999"
      }
      if (viajes2$utam_origen[i] %in% excl){
        viajes2$utam_origen[i] <- "UTAM999"
      }
    }

On fait sortir la personne de l’UTAM d’origine à l’heure de départ et on la fait rentrer dans l’UTAM de destination à l’heure d’arrivée.

    #Les heures et durées sont en fraction de la journée. On les multiplie par 24 pour les convertir en heures.
    viajes2$h_arrivee <- 24*viajes2$p31_hora_llegada
    viajes2$h_depart <- 24*viajes2$hora_inicio_viaje

    #pas de temps de l'animation : défaut 15 minutes
    pas_de_temps <- 0.25

    #On arrondit les heures au quart d'heure précédent. Par exemple 8h27 --> 8h15
    viajes2$h_arrivee_round <- pas_de_temps*(viajes2$h_arrivee-viajes2$h_arrivee%%pas_de_temps)%/%pas_de_temps
    viajes2$h_depart_round <- pas_de_temps*(viajes2$h_depart-viajes2$h_depart%%pas_de_temps)%/%pas_de_temps

On filtre pour ne garder que les trajets qui font une boucle dans la journée afin d’éviter un décrochage lors de la lecture en boucle de l’animation. Pour chaque trajet d’un individu, on calcule la somme UTAM d’arrivée - UTAM de départ, puis on somme ces valeurs pour chaque individu. Les individus qui font une boucle sont repérés par la valeur 0. Les autres, que l’on souhaite exclure, ont une valeur strictement positive ou négative. Cette opération retire 5.7% des observations, ce que l’on considère acceptable.

    viajes2$utam_origen_id <- as.numeric(gsub("UTAM", "", viajes2$utam_origen))
    viajes2$utam_destino_id <- as.numeric(gsub("UTAM", "", viajes2$utam_destino))
    viajes2$utam_test <- viajes2$utam_destino_id - viajes2$utam_origen_id

    viajes_circulares_2 <- viajes2 %>%
      group_by(id_hogar, id_persona) %>%
      summarize(circular_2 = sum(utam_test))

    viajes2 <- viajes2 %>%
      left_join(viajes_circulares_2, by = c("id_hogar" = "id_hogar", "id_persona" = "id_persona"))

    viajes2 <- viajes2 %>%
      filter(circular_2 == 0)

On retire les trajets intra-UTAM qui sont considérés, à cette échelle, comme des trajets n’ayant pas eu lieu car ils ne contribuent pas à la déformation.

    viajes2 <- viajes2 %>% filter(viajes2$utam_origen != viajes2$utam_destino)

On prépare un tableau contenant le solde de flux par UTAM

    #solde de mouvements par utam pour une heure donnée pondérée par le facteur d'expansion
    sorties_utam <- viajes2 %>%
      group_by(h_depart_round,utam_origen) %>%
      summarize(sorties = sum(f_exp))

    entrees_utam <- viajes2 %>%
      group_by(h_arrivee_round,utam_destino) %>%
      summarize(entrees = sum(f_exp))

    #création d'un tableau avec le solde de flux par UTAM
    flux <- data.frame(EMU2019$UTAM)
    names(flux)[names(flux) == 'EMU2019.UTAM'] <- 'UTAM'

    flux <- flux %>%
      full_join(sorties_utam, by = c("UTAM" = "utam_origen"))

    flux <- flux %>%  
      full_join(entrees_utam, by = c("UTAM" = "utam_destino",  "h_depart_round" = "h_arrivee_round"))

    flux[is.na(flux)] <- 0

    flux$solde <- flux$entrees - flux$sorties
    names(flux)[names(flux) == 'h_depart_round'] <- 'h_round'

    #tri chronologique par heure, puis par UTAM
    flux <- arrange(flux, h_round, UTAM)

    #conversion en matrice de flux
    f <- pivot_wider(flux, id_cols = UTAM, names_from = h_round, values_from = solde)
    f[is.na(f)] <- 0

Pour choisir l’heure de référence (où chacun est chez soi), on trace un histogramme de la répartition des heures de déplacement. On se rend compte que le creux de déplacement est maximal à 3h00. On va affecter la variable NUMPERSTOT, qui donne la population de nuit issue du recensement de 2018, à ce créneau.

    #entrées
    ggplot(flux, aes(x = h_round, y = entrees)) +
      stat_summary(geom = "bar", position = "dodge", fun = sum)

    #sorties 
    ggplot(flux, aes(x = h_round, y = sorties)) +
      stat_summary(geom = "bar", position = "dodge", fun = sum)

    #Création d'un tableau avec le stock de population par UTAM, heure par heure

    stock <- f %>%
      inner_join(EMU2019, by = 'UTAM', copy = TRUE) %>% #jointure avec EMU_2019
      relocate('NUMPERSTOT', .after = '2.75') #positionnement de NUMPERSTOT au niveau de la tranche 2h45 du matin

    #pour réordonner les colonnes en commençant à 3h00.
    stock <- stock[,1:98] 
    stock1 <- stock[,1]
    stock2 <- stock[,2:13]
    stock3 <- stock[,14:98]

    stock <-cbind(stock1, stock3, stock2)

    #remplissage de la matrice de stock
    for (i in 3:ncol(stock)){
      stock[,i]<-stock[,i]+stock[,i-1]
    }

    #On remplace les éventuelles valeurs négatives par des zéros. #Aucune UTAM concernée a priori
    for (i in 3:ncol(stock)){
      for (j in 1:nrow(stock)){
        if(stock[j,i]<0){
          stock[j,i] <- 0
        }
      }
    }

On transpose la table pour tracer des graphes de la population par UTAM au fil de la journée.

    stock2 <- stock
    rownames(stock2)<-stock2$UTAM #on renomme les lignes avec les codes des UTAM
    stock2 <- stock2[,3:98] #on supprime les deux premières colonnes
    tstock <- data.frame(t(stock2))

    tstock$hours <- rownames(tstock)
    tstock$rank <- c(1:96)

    #à titre d'exemple, on trace le graphe pour deux UTAM : EL LUCERO et CHAPINERO
    ggplot(tstock, aes(x = rank)) +                    
      geom_line(aes(y=UTAM67, col="EL LUCERO")) +  
      geom_line(aes(y=UTAM99, col="CHAPINERO")) + 
      scale_colour_manual("", breaks = c("EL LUCERO", "CHAPINERO"), values = c("red", "blue")) +
      scale_x_continuous(breaks=c(1+4*(0:23)), labels=c(3:23,0:2)) +
      scale_y_continuous(breaks=c(50000,100000), labels=c("50 000","100 000")) +
      labs(x="Time of day", y="Population") +
      theme(legend.position = "left")

Premier cartogramme

    #Récupération des valeurs de stock dans le fond de carte
    EMU2019_2 <- EMU2019[,c(2,55)] %>%
      left_join(stock, by = "UTAM")

    #on génère la liste des heures bien ordonnées
    heures <- c(0.25*(12:95),0.25*(0:11))
    #on génère des étiquettes à afficher sur les futures cartes (dans le titre) à partir de cette liste en les mettant au format anglophone
    heures2 <- heures
    for(i in 41:84){
      heures2[i] <- heures2[i]-12
    }
    for(i in 85:88){
      heures2[i] <- heures2[i]+12
    }
    
    heures_lab <- as.character(heures2)
    
    for(i in 1:length(heures_lab)){
      heures_lab[i] <- paste(as.character(heures2[i]%/%1), 
                         ".",as.character(60*heures2[i]%%1), sep = "")
      if(i%%4 == 1){heures_lab[i] <- str_replace(heures_lab[i], fixed(".0"), fixed(".00"))}
    }
    
    for(i in 1:length(heures_lab)){
      if(i %in% 37:84){
        heures_lab[i] <- paste(heures_lab[i], "pm", sep = " ")
      }
      else{
        heures_lab[i] <- paste(heures_lab[i], "am", sep = " ")
      }
    }
    
    #on renomme les colonnes des heures
    for (i in 1:ncol(EMU2019_2)){ 
      names(EMU2019_2)[names(EMU2019_2) == heures[i]] <- paste("stock",heures[i], sep = "_")
    }

Tests d’export de cartogrammes sur des heures individuelles, dans le cas présent à 3h00. On fixe une erreur absolue d’une valeur de 10 hectares, soit 100 000 m². On limite le calcul à 100 itérations.

    EMU2019_carto <- cartogramR(EMU2019_2, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 100, absrel = FALSE, abserror = 100000, L = 256))

    titre <- paste("Population per UTAM in Bogotá", heures_lab[1], sep = " - ")
    png(paste(titre,".png"), width = 960, height = 960)
    plot(EMU2019_carto, color = "black")
    title(main = titre, cex.main = 2)
    dev.off()

## png 
##   2

    include_graphics(paste(titre,".png"))

    #conversion en fichier sf
    EMU2019_carto_sf <- as.sf(EMU2019_carto)

Cartographie de la seconde variable : facteur multiplicatif de la population à chaque quart d’heure par rapport à la population nocturne

On va calculer les facteurs multiplicatifs à partir de la table de stock. On remplit le tableau colonne par colonne en divisant le stock pour chaque tranche horaire par le stock de référence NUMPERSTOT.

    variation <- stock

    for (i in 3:ncol(variation)){ 
      variation[,i]<-variation[,i]/variation[,2]
      #on renomme les colonnes pour les distinguer des variables de stock lors de la jointure
      names(variation)[names(variation) == heures[i]] <- paste("var",heures[i], sep = "_") 
    }
    names(variation)[names(variation) == heures[1]] <- paste("var",heures[1], sep = "_")
    names(variation)[names(variation) == heures[2]] <- paste("var",heures[2], sep = "_")

    #Jointure avec le fichier des géométries
    EMU2019_3 <- EMU2019_2 %>%
      left_join(variation, by = "UTAM")

Choix des seuils pour 6 classes

    bks <- c(0,0.7,0.9,1.1,2,5,50)

Création d’une palette de couleurs (double gradation harmonique)

    cols <- c("#303b98","#00beed","#fafaaa","#f2bf26","#e96831","#bb2b30")
    #facteur de transparence
    trans <- 50
    #application du facteur de transparence
    for (j in 1:6){
      cols[j] <- t_col(cols[j], percent = trans)
    }

Affichage de la carte chloroplèthe à 12h00, à titre d’exemple

      titre <- paste("Multiplication factor with respect to night population", heures_lab[37], sep = " - ")
      png(paste(titre, ".png"), width = 960, height = 960)
      plot(st_geometry(EMU2019_3))
      choroLayer(x = EMU2019_3, var = paste("var", heures[37], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1.5, legend.title.cex = 2)
      title(main = paste("Multiplication factor with respect to night population", heures_lab[37], sep = "\n"), cex.main = 2)
      dev.off()

## png 
##   2

      include_graphics(paste(titre,".png"))

Lissage spatial

Étapes préliminaires (calcul de l’autocorrélation spatiale et du variogramme)

    #Pour convertir la couche des UTAM au format spatial object (format requis par le package geoR)

    UTAM <- as(EMU2019_3, "Spatial")

    #Pour calculer les centroides des UTAM (indispensable pour le calcul sur les plus proches voisins)
    
    UTAMCentroids <- gCentroid(UTAM,byid=TRUE)

    #Pour récupérer les données initiales des UTAM sur les centroides

    UTAMCentroids <- SpatialPointsDataFrame(UTAMCentroids, UTAM@data)

    #Calcul sur les plus proches voisins

    listPPV <- knearneigh(UTAMCentroids@coords, k = 1) # pour connaître le plus proche voisin de chaque UTAM
    PPV <- knn2nb(listPPV, row.names = UTAM$UTAM) # pour convertir l'objet knn en objet nb
    distPPV <- nbdists(PPV, UTAMCentroids@coords) # pour connaître la distance entre plus proches voisins
    print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))

##       Min. 1st Qu.   Median     Mean  3rd Qu.     Max.
## 1 553.5302 1365.71 1587.616 2181.223 1854.208 10899.47

    hist(unlist(distPPV), breaks = 20,
         main = "Distance to the closest neighbour",
         col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")

La plupart des UTAM ont au moins un voisin dans un rayon de 1500m.

    #pour convertir les UTAM en objects nb
    nbUTAM <- poly2nb(pl = UTAM,
                     row.names = UTAM$UTAM,
                     snap = 50,
                     queen = TRUE)

    # pour identifier les UTAM sans voisins topologiques 
    summary(nbUTAM)

## Neighbour list object:
## Number of regions: 132 
## Number of nonzero links: 632 
## Percentage nonzero weights: 3.627181 
## Average number of links: 4.787879 
## 16 regions with no links:
## UTAM89 UTAM563 UTAM540 UTAM580 UTAM640 UTAM650 UTAM660 UTAM670 UTAM680
## UTAM690 UTAM700 UTAM600 UTAM630 UTAM620 UTAM610 UTAM590
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8  9 10 
## 16  3  4 12 19 19 25 16 11  6  1 
## 3 least connected regions:
## UTAM68 UTAM500 UTAM520 with 1 link
## 1 most connected region:
## UTAM112 with 10 links

    # création d'une liste des 16 UTAM sans voisins fournis par l'instruction précédente
    UTAM_isolees <- c("UTAM89",
              "UTAM563",
              "UTAM540",
              "UTAM580",
              "UTAM640",
              "UTAM650",
              "UTAM660",
              "UTAM670",
              "UTAM680",
              "UTAM690",
              "UTAM700",
              "UTAM600",
              "UTAM630",
              "UTAM620",
              "UTAM610",
              "UTAM590")

    #suppression des "UTAM_isolees" sans quoi le calcul de l'indice de Moran n'est pas possible
    UTAM <- UTAM[which(!UTAM$UTAM %in% UTAM_isolees),]
    EMU2019_3 <- EMU2019_3[which(!EMU2019_3$UTAM %in% UTAM_isolees),]
    nbUTAM <- poly2nb(pl = UTAM,
                     row.names = UTAM$UTAM,
                     snap = 50,
                     queen = TRUE)

    #Calcul du test de Moran pour chaque heure
    moran <- c(0*(1:96))
    for(i in 1:96){
      m <- moran.test(UTAM@data[,100+i], listw = nb2listw(nbUTAM)) #dans la table UTAM, les colonnes contenant les facteurs multiplicatifs sont les 101 à 196
      moran[i] <- as.numeric(m$estimate[1])
    }
    #affichage
    df <- data.frame(cbind(moran,heures))
    ggplot(df) +  
      theme_classic() +
      geom_point(aes(x = heures, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5)  +
      labs(x="Time of day", y="Moran") +
      theme(legend.position = "left", 
            panel.border = element_rect(colour = "black", fill=NA, size=0.5))

### Variogramme

    #pour calculer les pseudo-centroides des UTAM sans les UTAM isolées (indispensable pour le calcul du semi-variogramme)
    UTAMCentroids <- gPointOnSurface(UTAM,byid=TRUE)

    #pour récupérer les données initiales des communes sur les centroides
    UTAMCentroids <- SpatialPointsDataFrame(UTAMCentroids, UTAM@data)

    #pour convertir le SpatialPointsDataFrame en objet geodata
    UTAMCentroids.geodata <- as.geodata(UTAMCentroids, data.col = "var_12")

    #pour calculer le semi-variogramme empirique
    vario.ex<- variog(UTAMCentroids.geodata, bin.cloud=TRUE, option = "bin")

## variog: computing omnidirectional variogram

    plot(vario.ex, main = "Semivariogram of the multiplication factor at 12.00 pm in function of the distance", cex.main = 1)
    lines(vario.ex, type ="l", lty = 2, col="red")

### Carte lissée à 12h00

    #pour définir le contour de Bogotá comme emprise pour le lissage
    Emprise <- as.owin(gBuffer(UTAM, width=0))

    #pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
    UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$var_12)

    #pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha) --> calcul long
    cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))

    #Conversion de la surface lissée au format raster
    cartelissee.raster <- raster(cartelissee)
    crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster

    #Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
    par(mar = c(0, 0, 0, 0))

    #Affichage de toutes les UTAM en arrière plan en centrant la carte sur Bogotá.
    plot(st_geometry(EMU2019_3), col = NA, border = "black", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))

    #Affichage de la surface lissée
    plot(cartelissee.raster, breaks = bks, col=cols, add = T, legend=F)

    #Affichage des limites des UTAM
    plot(st_geometry(EMU2019_3), border = "black", lwd = 0.05, lty=3, add = T)

    legendChoro(
      pos = "bottomleft",
      title.txt = "Multiplication factor",
      breaks = bks, 
      nodata = FALSE,
      values.rnd = 2,
      col = cols,
      cex = 1.2,
      values.cex = 1,
      title.cex = 1
    )

    title("Population multiplication factor at 12.00 pm \n Smoothing radius 1000m", cex.main = 1, line = -2)

Le choix du rayon de lissage est important. Après plusieurs essais, on retient un rayon de 1000m qui est un bon compromis entre la prise en compte de suffisamment de voisins pour une UTAM donnée et la représentation la plus fidèle de la dynamique des différentes UTAM. Ci-dessous, la comparaison entre un rayon de lissage de 1000m et un rayon de 1500m.

##Isolement de l’UTAM “El Rincon de Suba” pour l’afficher en tant que légende sur le côté de la carte et représenter l’échelle de la déformation

  EMU2019_carto <- cartogramR(EMU2019_2, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 100, absrel = FALSE, abserror = 100000, L = 256))
  EMU2019_carto_sf <- as.sf(EMU2019_carto)
  RINCON <- EMU2019_carto_sf[which(EMU2019_carto_sf$UTAM == "UTAM28"),]
  st_geometry(RINCON)[[1]] = st_geometry(RINCON)[[1]] + st_point(c(-25000,-10000))
  RINCON$legend <- paste(RINCON$NUMPERSTOT, "people", sep = "\n")

Déformation initiale de la carte pour le cartogramme seul

Pour ce faire, on réalise une interpolation entre la carte de Bogotá non déformée et la carte déformée selon la population de nuit à 3h00.

  #Déformation du fond de carte initial pour le cartogramme seul
    dir.create("Animated_Cartogram")

    #Calcul des surfaces
    EMU2019_area <- EMU2019[which(!EMU2019$UTAM %in% UTAM_isolees),]

    EMU2019_area$Area <- as.numeric(st_area(st_geometry(EMU2019_area)))

    #calcul de la population totale
    Population <- sum(EMU2019_area$NUMPERSTOT)

    #calcul de la superficie totale
    Superficie <- sum(EMU2019_area$Area)

    #nombre d'images intermédiaires pour la déformation
    imgs <- 8

    #Ventilation de la population totale proportionnellement à la superficie des UTAM pour l'image 0
    EMU2019_area$NUMPERSTOT_0 <- EMU2019_area$Area*Population/Superficie

    #Champ que l'on va faire varier
    EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT

    #création des images intermédiaires par une interpolation linéaire (plus lisse qu'en utilisant le TCAM, qui donnerait une croissance trop rapide sur la fin et un arrêt brutal)
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    par(mar = c(1, 1, 1, 1))

    #affichage de la carte non déformée
    png(paste0(getwd(),"/Animated_Cartogram/", "deformation 0.png"), width = 960, height = 960)
    plot(st_geometry(EMU2019_area), border = "black", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
    dev.off()

## png 
##   2

    for (i in 1:imgs){
      EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT*i/imgs + EMU2019_area$NUMPERSTOT_0*(imgs-i)/imgs
      carto <- cartogramR(EMU2019_area, count = "NUMPERSTOT_i", options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      png(paste0(getwd(),"/Animated_Cartogram/", paste0("deformation ", i,".png")), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = NA, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      carto_sf <- as.sf(carto)
      plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
      dev.off()
    }

Déformation initiale de la carte pour le cartogramme avec seconde variable (lissée ou non, discrète ou continue)

  #Déformation du fond de carte initial pour le cartogramme seconde variable. 
    dir.create("Choropleth_animated_cartogram")

    #Calcul des surfaces
    EMU2019_area <- EMU2019[which(!EMU2019$UTAM %in% UTAM_isolees),]

    EMU2019_area$Area <- as.numeric(st_area(st_geometry(EMU2019_area)))

    #calcul de la population totale
    Population <- sum(EMU2019_area$NUMPERSTOT)

    #calcul de la superficie totale
    Superficie <- sum(EMU2019_area$Area)

    #nombre d'images intermédiaires pour la déformation
    imgs <- 8

    #Ventilation de la population totale proportionnellement à la superficie des UTAM pour l'image 0
    EMU2019_area$NUMPERSTOT_0 <- EMU2019_area$Area*Population/Superficie

    #Champ que l'on va faire varier
    EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT

    #création des images intermédiaires par une interpolation linéaire (plus lisse qu'en utilisant le TCAM, qui donnerait une croissance trop rapide sur la fin et un arrêt brutal)
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    par(mar = c(1, 1, 1, 1))

    #affichage de la carte non déformée
    png(paste0(getwd(),"/Choropleth_animated_cartogram/", "deformation 0.png"), width = 960, height = 960)
    plot(st_geometry(EMU2019_area), col = NA, border = "black", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))

    dev.off()

## png 
##   2

    for (i in 1:imgs){
      EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT*i/imgs + EMU2019_area$NUMPERSTOT_0*(imgs-i)/imgs
      carto <- cartogramR(EMU2019_area, count = "NUMPERSTOT_i", options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      png(paste0(getwd(),"/Choropleth_animated_cartogram/", paste0("deformation ", i,".png")), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      carto_sf <- as.sf(carto)
      plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
      dev.off()
    }
    
    # duplicate first initial images
    dir.create("Smoothed_binned_colored_animated_cartogram")
    initial_files <- list.files(paste0(getwd(),"/Choropleth_animated_cartogram/"))
    file.copy(from = paste0(getwd(),"/Choropleth_animated_cartogram/", initial_files),
          to = paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram/", initial_files))

## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

    dir.create("Smoothed_continuous_colored_animated_cartogram")
    file.copy(from = paste0(getwd(),"/Choropleth_animated_cartogram/", initial_files),
          to = paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram/", initial_files))

## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

Fondu enchaîné pour l’apparition de la légende pour la carte affichant le cartogramme seul

    setwd(paste0(getwd(),"/Animated_Cartogram"))
    carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
  
  #niveau de transparence (%)
  
  for (i in c(90, 80, 70, 60, 50)) {
    trans <- i
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    par(mar = c(1, 1, 1, 1))
    
    # Affichage de toutes les UTAM en arrière-plan
    titre <- paste("transition", 100 - trans, "percent", sep = " - ")
    png(paste(titre,".png"), width = 960, height = 960)
    plot(st_geometry(EMU2019_3), col = NA, border = "white", bg = "white", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
    
    par(col.main = paste("grey", 1.8*i-90))
    title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2) 
    
    # conversion de l'objet cartogram en objet sf
    carto_sf <- as.sf(carto)
    plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
    
    #Affichage de l'UTAM EL RINCON dans un coin
    plot(st_geometry(RINCON), col = NA, border = paste("grey", 1.8*i-90), add = TRUE)
    labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0, col = paste("grey", 1.8*i-90))
    
    #Affichage de la carte non déformée
    mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
    split.screen(figs = mat, erase=FALSE)
    plot(st_geometry(EMU2019_3), border = paste("grey", 0.6*i+30), add = TRUE)
    mtext("Initial undistorted map", side = 1, line = 3, col = paste("grey", 1.8*i-90), outer = FALSE)
    close.screen(all = TRUE)
    dev.off()
  }

Fondu enchaîné pour l’apparition de la légende et de la seconde variable pour la carte affichant la seconde variable non lissée

    setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))
    
    Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"
    
    carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
    
    EMU2019_carto_sf <- as.sf(carto)
    
    for (i in c(90,80,70,60,50)) { 
      trans <- i
      #Affichage du fond
      par(mar = c(1,1,1,1))
      titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      choroLayer(x = EMU2019_carto_sf, var = paste("var", heures[1], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1, legend.title.cex = 1.5, add = TRUE)
      
      title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2) 
      
      #Affichage de l'UTAM EL RINCON dans un coin
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
      
      #Affichage de la carte non déformée
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(n = 1)
      
      background <- t_col("grey90",200-2*i)
      mat <- matrix(c(0,1.0,0,1), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), col = NA, border = background, bg = background, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      plot(st_geometry(EMU2019_carto_sf), border = "black", lwd = 1, add = TRUE)
      close.screen(all = TRUE)
      dev.off()
    }

Fondu enchaîné pour l’apparition de la légende et de la seconde variable pour la carte affichant le cartogramme et la surface lissée discrétisée

    setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
    
    Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"
    
    carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))

    # pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
    UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,196])
    
    # pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha) --> calcul long
    cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
    
    # Conversion de la surface lissée au format raster
    cartelissee.raster <- raster(cartelissee)
    crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
    
    # reclassification de la surface lissée préalablement calculée
    cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
    
    #vectorisation de la surface reclassée
    cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
    
    # pour trier les enregistrements de la colonne layer
    cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
    
    # création d'une liste de classes présentes à l'heure donnée
    liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
    
    # Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
    cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
    
    # Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
    cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
    
    # pour généraliser les contours de cartelissee.vecteur.drapee (on garde 5% des sommets initiaux)
    cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.05, keep_shapes = TRUE)
    
    #niveau de transparence (%)
    
    for (i in c(90, 80, 70, 60, 50)) {
      trans <- i
      cols_2 <- cols
      for (j in 1:10){
        cols_2[j] <- t_col(cols[j], percent = trans)
      }
      
      Emprise <- as.owin(gBuffer(UTAM, width=0))
      par(mar = c(1, 1, 1, 1))
      
      # création d'un tableau de correspondances avec pour chaque classe son code couleur
      Correspondance_palette <- as.data.frame(cbind(c(1:6), cols_2))
      names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"
      
      # création d'une palette de couleurs correspondant uniquement aux classes présentes
      couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
      
      # Affichage de toutes les UTAM en arrière-plan
      titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      
      # Affichage de la surface lissée déformée
      typoLayer(
        x = cartelissee.vecteur.drapee,
        var="layer",
        col = as.character(couleurs_classes_presentes$cols),
        lwd = 0.1,
        border = as.character(couleurs_classes_presentes$cols),
        legend.pos = "n",
        add=T)
      
      legendChoro(
        pos = "bottomleft",
        title.txt = "Multiplication factor",
        breaks = c("","0.72","0.9","1.1","2","5",""),
        nodata = FALSE,
        values.rnd = 2,
        col = cols,
        cex = 1.2,
        values.cex = 1,
        title.cex = 1.5
      )
      
      par(col.main = paste("grey", 1.8*i-90))
      title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2) 
      
      # conversion de l'objet cartogram en objet sf
      carto_sf <- as.sf(carto)
      
      # Affichage
      plot(st_geometry(carto_sf), border = paste("grey", 200-2*i), lwd = 1, add = TRUE)
      
      #Affichage des contours des localités
      Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
      Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
      plot(st_geometry(Localidades_carto_2), col = NA, border = "black", lwd = 3.25 - i/40, add = TRUE)
      
      
      #Affichage de l'UTAM EL RINCON en légende
      plot(st_geometry(RINCON), col = NA, border = paste("grey", 1.8*i-90), add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0, col = paste("grey", 1.8*i-90))
      
      #Affichage de la carte non déformée dans un carton en bas à droite
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = paste("grey", 0.6*i+30), add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, col = paste("grey", 1.8*i-90), outer = FALSE)
      close.screen(all = TRUE)
      dev.off()
    }

Fondu enchaîné pour l’apparition de la légende et de la seconde variable pour la carte affichant le cartogramme et la surface lissée continue

    setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
    # paramétrage de la discrétisation (définition des seuils ou bornes de classes)
    bks_pseudo_continuous <- c(seq(from=0, to=3, by = 0.025), seq(from=3.5, to=50, by = 0.5))
    
    n0<-length(seq(from=0, to=0.6, by = 0.025))
    n1<-length(seq(from=0.625, to=0.8, by = 0.025))
    n2<-length(seq(from=0.825, to=0.975, by = 0.025))
    n3<-length(seq(from=1, to=1.5, by = 0.025))
    n4<-length(seq(from=1.525, to=3, by = 0.025))
    n5<-length(seq(from=3.5, to=10, by = 0.5))
    n6<-length(seq(from=10.5, to=50, by = 0.5))
    
    
    #palette de couleur (double gradation harmonique) et facteur de transparence
    cols_continuous <- c("#1d235c", "#303b98","#00beed","#fafaaa","#f2bf26","#e96831","#bb2b30","#7f0000")
    trans <- 50
    
    pal0 <- colorRampPalette(colors = cols_continuous[1:2], interpolate = "linear")
    palette0 <- pal0(n0)
    palette0 <- t_col_palette(palette0, trans)
    
    pal1 <- colorRampPalette(colors = cols_continuous[2:3], interpolate = "linear")
    palette1 <- pal1(n1)
    palette1 <- t_col_palette(palette1, trans)
    
    pal2 <- colorRampPalette(colors = cols_continuous[3:4], interpolate = "linear")
    palette2 <- pal2(n2)
    palette2 <- t_col_palette(palette2, trans)
    
    pal3 <- colorRampPalette(colors = cols_continuous[4:5], interpolate = "linear")
    palette3 <- pal3(n3)
    palette3 <- t_col_palette(palette3, trans)
    
    pal4 <- colorRampPalette(colors = cols_continuous[5:6], interpolate = "linear")
    palette4 <- pal4(n4)
    palette4 <- t_col_palette(palette4, trans)
    
    pal5 <- colorRampPalette(colors = cols_continuous[6:7], interpolate = "linear")
    palette5 <- pal5(n5)
    palette5 <- t_col_palette(palette5, trans)
    
    pal6 <- colorRampPalette(colors = cols_continuous[7:8], interpolate = "linear")
    palette6 <- pal6(n6)
    palette6 <- t_col_palette(palette6, trans)
    
    palette <- c(palette0, palette1, palette2, palette3, palette4, palette5, palette6)
    
    #preparation de la légende continue
    leg1 <- pal1(10)
    leg1 <- t_col_palette(leg1, trans)
    leg2 <- pal2(10)
    leg2 <- t_col_palette(leg2, trans)
    leg3 <- pal3(10)
    leg3 <- t_col_palette(leg3, trans)
    leg4 <- pal4(10)
    leg4 <- t_col_palette(leg4, trans)
    leg5 <- pal5(10)
    leg5 <- t_col_palette(leg5, trans)
    leg <- cbind(leg1, leg2, leg3, leg4, leg5)
    
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    
    Correspondance_palette <- as.data.frame(cbind(c(1:(n0+n1+n2+n3+n4+n5+n6)), palette))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V2"] <- "palette"
    
    carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
    
    #pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
    UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,196])
    
    #pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
    cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
    
    #Conversion de la surface lissée au format raster
    cartelissee.raster <- raster(cartelissee)
    crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
    
    # pour découper le raster sur l'emprise des UTAM
    cartelissee.raster <- mask(cartelissee.raster, UTAM)
    
    #reclassification de la surface lissée préalablement calculée
    cartelissee.reclass <- cut(cartelissee.raster, breaks = bks_pseudo_continuous)
    
    #vectorisation de la surface reclassée
    cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
    
    #pour trier les enregistrements de la colonne layer
    cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
    
    #création d'une liste de classes présentes à l'heure donnée
    liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
    
    l0 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer <= n0]
    l1 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+1):(n0+n1)]
    l2 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+1):(n0+n1+n2)]
    l3 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+1):(n0+n1+n2+n3)]
    l4 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+1):(n0+n1+n2+n3+n4)]
    l5 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+n4+1):(n0+n1+n2+n3+n4+n5)]
    l6 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer > (n0+n1+n2+n3+n4+n5)]
    
    #création d'une palette de couleurs correspondant uniquement aux classes présentes
    couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
    
    # Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
    cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
    
    #Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
    cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
    
    #pour généraliser les contours de cartelissee.vecteur.drapee (on garde 10% des sommets initiaux)
    cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.1, keep_shapes = TRUE)
    
    #niveau de transparence (%)
    
    for (i in c(50,60,70,80,90)) { 
      trans <- i
      #Affichage du fond
      par(mar = c(1,1,1,1))
      titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      
      #Affichage de la surface lissée déformée (valeurs jusqu'à 0.6)
      if(length(l0)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l0, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l0]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeurs de 0.6 à 0.8)
      if(length(l1)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l1, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l1]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 0.8 à 1)
      if(length(l2)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l2, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l2]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 1 à 1.5)
      if(length(l3)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l3, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l3]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 1.5 à 3)
      if(length(l4)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l4, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l4]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 3 à 10)
      if(length(l5)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l5, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l5]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur à partir de 10)
      if(length(l6)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l6, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l6]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #seulement pour l'échelle de la légende, sans les couleurs
      legendChoro(
        pos = "bottomleft",
        title.txt = "Multiplication factor",
        breaks = c("","0.72","0.9","1.1","2","5",""),
        nodata = FALSE,
        values.rnd = 2,
        col = c("grey90","grey90","grey90","grey90","grey90","grey90"),
        cex = 1.2,
        values.cex = 1,
        title.cex = 1.5,
        border = "grey90"
      )
      
      title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2) 
      
      # conversion de l'objet cartogram en objet sf
      carto_sf <- as.sf(carto)
      
      #Affichage des contours des localités
      Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
      Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
      
      #Affichage de l'UTAM EL RINCON dans un coin
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
      
      #affichage de la légende continue
      mat <- matrix(c(0.013,0.14,0.013,0.297), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      
      image(1, 1:length(leg), t(as.matrix(1:length(leg))), col = leg, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n",bty = "n")
      close.screen(n = 1)
      
      #Affichage de la carte non déformée
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(n = 1)
      
      background <- t_col("grey90",200-2*i)
      mat <- matrix(c(0,1.0,0,1), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), col = NA, border = background, bg = background, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      plot(st_geometry(carto_sf), border = paste("grey", 200-2*i), lwd = 1, add = TRUE)
      plot(st_geometry(Localidades_carto_2), col = NA, border = "black", lwd = 3.25 - i/40, add = TRUE)
      close.screen(all = TRUE)
      dev.off()
    }

Cartogramme seul au pas de temps d’un quart d’heure (création de la série d’images)

    setwd(paste0(getwd(),"/Animated_Cartogram"))
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    
    for (i in 1:96){
      par(mar = c(1, 1, 1, 1))
      carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      
      # Affichage du fond
      plot(st_geometry(EMU2019_3), col = NA, border = "white", bg = "white", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      
      title(main = titre, cex.main = 2, line = -2) 
      
      # conversion de l'objet cartogram en objet sf
      carto_sf <- as.sf(carto)
      plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
      
      #Affichage de l'UTAM EL RINCON en légende
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
      
      #Affichage de la carte non déformée dans un carton en bas à droite
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(all = TRUE)
      dev.off()
    }

Cartogramme et seconde variable non lissée (création de la série d’images)

    setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))

    Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"

    for (i in 1:96){
      carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      EMU2019_carto_sf <- as.sf(carto)
      par(mar = c(1, 1, 1, 1))
      titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      choroLayer(x = EMU2019_carto_sf, var = paste("var", heures[i], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1, legend.title.cex = 1.5, add = TRUE)
      title(main = titre, cex.main = 2, line = -2) 
      #Affichage de l'UTAM EL RINCON dans un coin
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
      #Affichage de la carte non déformée
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(all = TRUE)
      dev.off()
    }

Cartogramme et surface lissée discrétisée au pas de temps d’un quart d’heure (création de la série d’images)

    setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"

    for (i in 1:96){
      par(mar = c(1, 1, 1, 1))
      carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)

      #pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
      UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,100+i])

      #pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
      cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))

      #Conversion de la surface lissée au format raster
      cartelissee.raster <- raster(cartelissee)
      crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster

      #reclassification de la surface lissée préalablement calculée
      cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)

      #vectorisation de la surface reclassée
      cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE)

      #pour trier les enregistrements de la colonne layer
      cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]

      #création d'une liste de classes présentes à l'heure donnée
      liste_classes_presentes <- as.character(cartelissee.vecteur$layer)

      #création d'une palette de couleurs correspondant uniquement aux classes présentes
      couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
      
      # Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
      cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)

      #Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
      cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)

      #pour généraliser les contours de cartelissee.vecteur.drapee (on garde 5% des sommets initiaux)
      cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.05, keep_shapes = TRUE)

      #Affichage du fond
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))

      #Affichage de la surface lissée déformée
      typoLayer(
        x = cartelissee.vecteur.drapee,
        var="layer",
        col = as.character(couleurs_classes_presentes$cols),
        lwd = 0.1,
        border = NA,
        legend.pos = "n",
        add=T)

      legendChoro(
        pos = "bottomleft",
        title.txt = "Multiplication factor",
        breaks = c("","0.72","0.9","1.1","2","5",""),
        nodata = FALSE,
        values.rnd = 2,
        col = cols,
        cex = 1.2,
        values.cex = 1,
        title.cex = 1.5
      )

      title(main = titre, cex.main = 2, line = -2) 

      #conversion de l'objet cartogram en objet sf
      carto_sf <- as.sf(carto)
      
      #Affichage des limites des UTAM déformées
      plot(st_geometry(carto_sf), border = "white", lwd = 1, add = TRUE)
      
      #Affichage des zones d'étude MODURAL (optionnel)
      #plot(st_geometry(carto_sf_modural), border = "black", lwd = 3, add = TRUE)
      #plot(st_geometry(carto_sf_modural), border = "white", lwd = 1, add = TRUE)

      #Affichage des contours des localités
      Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
      Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
      plot(st_geometry(Localidades_carto_2), col = NA, lwd = 2, add = TRUE)

      #Affichage de l'UTAM EL RINCON en légende
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)

      #Affichage de la carte non déformée dans un carton en bas à droite
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(all = TRUE)
      dev.off()
    }

Cartogramme et surface lissée continue au pas de temps d’un quart d’heure (création de la série d’images)

    setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
    
    Emprise <- as.owin(gBuffer(UTAM, width=0))
    
    Correspondance_palette <- as.data.frame(cbind(c(1:(n0+n1+n2+n3+n4+n5+n6)), palette))
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V1"] <- "Id_Classe"
    names(Correspondance_palette)[names(Correspondance_palette) ==  "V2"] <- "palette"

    for (i in 1:96){
      par(mar = c(1, 1, 1, 1))
      carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
      titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
      png(paste(titre,".png"), width = 960, height = 960)
      
      #pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
      UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,100+i])
      
      #pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
      cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
      
      #Conversion de la surface lissée au format raster
      cartelissee.raster <- raster(cartelissee)
      crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
      
      # pour découper le raster sur l'emprise des UTAM
      cartelissee.raster <- mask(cartelissee.raster, UTAM)
      
      #reclassification de la surface lissée préalablement calculée
      cartelissee.reclass <- cut(cartelissee.raster, breaks = bks_pseudo_continuous)
      
      #vectorisation de la surface reclassée
      cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
      
      #pour trier les enregistrements de la colonne layer
      cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
      
      #création d'une liste de classes présentes à l'heure donnée
      liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
      
      l0 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer <= n0]
      l1 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+1):(n0+n1)]
      l2 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+1):(n0+n1+n2)]
      l3 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+1):(n0+n1+n2+n3)]
      l4 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+1):(n0+n1+n2+n3+n4)]
      l5 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+n4+1):(n0+n1+n2+n3+n4+n5)]
      l6 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer > (n0+n1+n2+n3+n4+n5)]
      
      #création d'une palette de couleurs correspondant uniquement aux classes présentes
      couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
      
      # Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
      cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
      
      #Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
      cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
      
      #pour généraliser les contours de cartelissee.vecteur.drapee (on garde 10% des sommets initiaux)
      cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.1, keep_shapes = TRUE)
      
      #Affichage du fond
      plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
      
      #Affichage de la surface lissée déformée (valeurs jusqu'à 0.6)
      if(length(l0)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l0, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l0]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeurs de 0.6 à 0.8)
      if(length(l1)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l1, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l1]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 0.8 à 1)
      if(length(l2)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l2, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l2]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 1 à 1.5)
      if(length(l3)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l3, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l3]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 1.5 à 3)
      if(length(l4)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l4, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l4]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur de 3 à 10)
      if(length(l5)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l5, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l5]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #Affichage de la surface lissée déformée (valeur à partir de 10)
      if(length(l6)>0){
        typoLayer(
          x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l6, ],
          var="layer",
          col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l6]),
          lwd = 0.1,
          border = NA,
          legend.pos = "n",
          add=T)
      }
      
      #seulement pour l'échelle de la légende, sans les couleurs
      legendChoro(
        pos = "bottomleft",
        title.txt = "Multiplication factor",
        breaks = c("","0.72","0.9","1.1","2","5",""),
        nodata = FALSE,
        values.rnd = 2,
        col = c("grey90","grey90","grey90","grey90","grey90","grey90"),
        cex = 1.2,
        values.cex = 1,
        title.cex = 1.5,
        border = "grey90"
      )
      
      title(main = titre, cex.main = 2, line = -2) 
      
      # conversion de l'objet cartogram en objet sf
      carto_sf <- as.sf(carto)
      # Affichage
      plot(st_geometry(carto_sf), border = "white", lwd = 1, add = TRUE)
      
      #Affichage des contours des localités
      Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
      Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
      plot(st_geometry(Localidades_carto_2), col = NA, lwd = 2, add = TRUE)
      
      #Affichage de l'UTAM EL RINCON dans un coin
      plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
      labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
      
      #affichage de la légende continue
      mat <- matrix(c(0.013,0.14,0.013,0.297), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      
      image(1, 1:length(leg), t(as.matrix(1:length(leg))), col = leg, xlab = "", ylab = "", xaxt = "n",
            yaxt = "n",bty = "n")
      close.screen(n = 1)
      
      #Affichage de la carte non déformée
      mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
      split.screen(figs = mat, erase=FALSE)
      plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
      mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
      close.screen(all = TRUE)
      dev.off()
    }

Création d’un gif animé à partir de la série d’images produites (cartogramme seul)

    dir.create("Animated_Cartogram_Export")
    
    #chargement des informations des fichiers à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Animated_Cartogram"), full.names = TRUE))
    
    #tri par date de création
    details = details[with(details, order(as.POSIXct(ctime))), ] 
    
    list <- rep("0",110)
    for (i in 1:10){
      list[i] <- paste("00",i, sep = "")
    }
    
    for (i in 10:110){
      if (i>=10 & i<100){
        list[i] <- paste("0",i, sep = "")
      } else {
        list[i] <- paste(i)
      }
    }
    
    file.rename(rownames(details), paste0(getwd(),"/Animated_Cartogram/anim", list[1:110],".png"))

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE

    # rechargement des informations des fichiers renommés à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Animated_Cartogram"), full.names = TRUE))
    
    details <- details[15:110,]

    imgs = rownames(details)

    img_list <- lapply(imgs, image_read)

    img_joined <- image_join(img_list)
    img_animated <- image_animate(img_joined, fps = 20)

    #export en gif
    image_write(image = img_animated, path = "Animated_Cartogram_Export/animation.gif")
    include_graphics(paste0(getwd(),"/Animated_Cartogram_Export/animation.gif"))

Création d’une vidéo mp4 à partir de la série d’images produites (cartogramme seul)

    #export en mp4
    rootwd <- getwd()
    setwd(paste0(getwd(),"/Animated_Cartogram"))
    av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")

## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Animated_Cartogram\\animation.mp4"

    file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Animated_Cartogram_Export"))

Création d’un gif animé à partir de la série d’images produites (cartogramme avec seconde variable non lissée)

    dir.create("Choropleth_animated_cartogram_Export")
    
    #chargement des informations des fichiers à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Choropleth_animated_cartogram"), full.names = TRUE))
    
    #tri par date de création
    details = details[with(details, order(as.POSIXct(ctime))), ] 
    
    list <- rep("0",110)
    for (i in 1:10){
      list[i] <- paste("00",i, sep = "")
    }
    
    for (i in 10:110){
      if (i>=10 & i<100){
        list[i] <- paste("0",i, sep = "")
      } else {
        list[i] <- paste(i)
      }
    }
    
    file.rename(rownames(details), paste0(getwd(),"/Choropleth_animated_cartogram/anim", list[1:110],".png"))

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE

    # rechargement des informations des fichiers renommés à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Choropleth_animated_cartogram"), full.names = TRUE))
    
    details <- details[15:110,]

    imgs = rownames(details)

    img_list <- lapply(imgs, image_read)

    img_joined <- image_join(img_list)
    img_animated <- image_animate(img_joined, fps = 20)

    #export en gif
    image_write(image = img_animated, path = "Choropleth_animated_cartogram_Export/animation.gif")
    include_graphics(paste0(getwd(),"/Choropleth_animated_cartogram_Export/animation.gif"))

Création d’une vidéo mp4 à partir de la série d’images produites (cartogramme avec seconde variable non lissée)

    #export en mp4
    rootwd <- getwd()
    setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))
    av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")

## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Choropleth_animated_cartogram\\animation.mp4"

    file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Choropleth_animated_cartogram_Export"))

Création d’un gif animé à partir de la série d’images produites (cartogramme et seconde variable lissée discrétisée)

    dir.create("Smoothed_binned_colored_animated_cartogram_Export")
    
    #chargement des informations des fichiers à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"), full.names = TRUE))
    
    #tri par date de création
    details = details[with(details, order(as.POSIXct(ctime))), ] 
    
    list <- rep("0",110)
    for (i in 1:10){
      list[i] <- paste("00",i, sep = "")
    }
    
    for (i in 10:110){
      if (i>=10 & i<100){
        list[i] <- paste("0",i, sep = "")
      } else {
        list[i] <- paste(i)
      }
    }
    
    file.rename(rownames(details), paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram/anim", list[1:110],".png"))

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE

    # rechargement des informations des fichiers renommés à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"), full.names = TRUE))
        
    details <- details[15:110,]

    imgs = rownames(details)

    img_list <- lapply(imgs, image_read)

    img_joined <- image_join(img_list)
    img_animated <- image_animate(img_joined, fps = 20)
        
    #export en gif
    image_write(image = img_animated, path = "Smoothed_binned_colored_animated_cartogram_Export/animation.gif")
    
    include_graphics(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram_Export/animation.gif"))

Création d’une vidéo mp4 à partir de la série d’images produites (cartogramme et seconde variable discrétisée lissée)

    #export en mp4
    rootwd <- getwd()
    setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
    av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")

## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Smoothed_binned_colored_animated_cartogram\\animation.mp4"

    file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Smoothed_binned_colored_animated_cartogram_Export"))

Création d’un gif animé à partir de la série d’images produites (cartogramme et seconde variable lissée continue)

    dir.create("Smoothed_continuous_colored_animated_cartogram_Export")
    
    #chargement des informations des fichiers à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"), full.names = TRUE))
    
    #tri par date de création
    details = details[with(details, order(as.POSIXct(ctime))), ] 
    
    list <- rep("0",110)
    for (i in 1:10){
      list[i] <- paste("00",i, sep = "")
    }
    
    for (i in 10:110){
      if (i>=10 & i<100){
        list[i] <- paste("0",i, sep = "")
      } else {
        list[i] <- paste(i)
      }
    }
    
    file.rename(rownames(details), paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram/anim", list[1:110],".png"))

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE

    # rechargement des informations des fichiers renommés à intégrer à l'animation
    details = file.info(list.files(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"), full.names = TRUE))
      
    details <- details[15:110,]

    imgs = rownames(details)

    img_list <- lapply(imgs, image_read)

    img_joined <- image_join(img_list)
    img_animated <- image_animate(img_joined, fps = 20)

    #export en gif
    image_write(image = img_animated, path = "Smoothed_continuous_colored_animated_cartogram_Export/animation.gif")
    include_graphics(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram_Export/animation.gif"))

Création d’une vidéo mp4 à partir de la série d’images produites (cartogramme et seconde variable lissée continue)

    #export en mp4
    rootwd <- getwd()
    setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
    av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")

## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Smoothed_continuous_colored_animated_cartogram\\animation.mp4"

    file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Smoothed_continuous_colored_animated_cartogram_Export"))